doseresNMA antidep/drnma.tab.comp.R

In htx-r/doseresNMA:

#!! not tested yet
require(R2jags)
require(coda)
require(dplyr)
# table of tau, %drop in tau2, g(reg coeff) and DIC
tab.comp <- function(mod1=drnma_RE,
                     mod2=drnma_rob,
                     mod3=drnma_year,
                     mod4=drnma_var,
                     mod5=drnma_class){
  # 5 models
  m1 <- round(mod1$BUGSoutput$summary,3)
  m2 <- round(mod2$BUGSoutput$summary,3)
  m3 <- round(mod3$BUGSoutput$summary,3)
  m4 <- round(mod4$BUGSoutput$summary,3)
  m5 <- round(mod5$BUGSoutput$summary,3)
  
  # tau (mean and  95%CrI)
  tau <- rbind(
           paste0(m1['tau.d','mean'],'[', m1['tau.d','2.5%'],',', m1['tau.d','97.5%'],']'),
           paste0(m2['tau.d','mean'],'[', m2['tau.d','2.5%'],',', m2['tau.d','97.5%'],']'),
           paste0(m3['tau.d','mean'],'[', m3['tau.d','2.5%'],',', m3['tau.d','97.5%'],']'),
           paste0(m4['tau.d','mean'],'[', m4['tau.d','2.5%'],',', m4['tau.d','97.5%'],']'),
           paste0(m5['tau.d','mean'],'[', m5['tau.d','2.5%'],',', m5['tau.d','97.5%'],']'))
  # drop in tau2
  tau.mean <- c(m1['tau.d','mean'],
                m2['tau.d','mean'],
                m3['tau.d','mean'],
                m4['tau.d','mean'],
                m5['tau.d','mean'])
  tau.mean1 <- m1['tau.d','mean']
  tau.mean <- m5['tau.d','mean']
  drop_tau2 <-  round(100*(tau.mean1^2-tau.mean^2)/tau.mean1^2)
  drop_tau2 <- round(100*(tau.mean[1]^2-tau.mean^2)/tau.mean[1]^2)
  
  # g - metareg coef (mean and 95%CrI)
  g <- rbind(
             paste0(m2['g','mean'],'[', m2['g','2.5%'],',', m2['g','97.5%'],']'),
             paste0(m3['g','mean'],'[', m3['g','2.5%'],',', m3['g','97.5%'],']'),
             paste0(m4['g','mean'],'[', m4['g','2.5%'],',', m4['g','97.5%'],']')
             )
  exp.g <- rbind(
    paste0(round(exp(m2['g','mean']),2),'[', round(exp(m2['g','2.5%']),2),',', round(exp(m2['g','97.5%']),2),']'),
    paste0(round(exp(m3['g','mean']),2),'[', round(exp(m3['g','2.5%']),2),',', round(exp(m3['g','97.5%']),2),']'),
    paste0(round(exp(m4['g','mean']*0.1),2),'[', round(exp(m4['g','2.5%']*0.1),2),',', round(exp(m4['g','97.5%']*0.1),2),']')
  )
  # DIC
  DIC1 <- dic.fun(mod1)
  DIC2 <- dic.fun(mod2)
  DIC3 <- dic.fun(mod3)
  DIC4 <- dic.fun(mod4)
  DIC5 <- dic.fun(mod5)
  
  DIC <- c(DIC1,
           DIC2,
           DIC3,
           DIC4,
           DIC5)
  # return: table
    mod.comp.tbl <- tibble(tau=tau,
                           drop_tau2= drop_tau2,
                           g=g,
                           DIC=DIC 
    )
    mod.comp.tbl
    
}




dic.fun <- function(x){
  jagssamples <- as.mcmc(x) # I NEED R2jags and coda packages to run that
  samples = do.call(rbind, jagssamples)%>% data.frame()
  dev <- samples %>% select(., starts_with("resdev"))
  totresdev <- samples$totresdev %>% mean()
  pmdev <- colMeans(dev)
  rhat <- samples %>% select(., starts_with("rhat"))
  r <- samples %>% select(., starts_with("r."))
  n <- samples %>% select(., starts_with("n"))
  rtilde <- rhat %>%
    colMeans() %>%
    matrix(nrow=1, ncol=ncol(rhat)) %>%
    data.frame() %>%
    slice(rep(1:n(), each = nrow(rhat)))
  
  pmdev_fitted <- 2*(r*log(r/rtilde)+(n-r)*log((n-r)/(n-rtilde)))[1,]
  leverage = pmdev-pmdev_fitted
  DIC = sum(leverage,na.rm = T) + totresdev
  DIC
}